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_■. ■ Molecular dynamics simulation is used to investigate the crystallization of a classi- 

J> . cal two-dimensional electron system, in which electrons interact with the Coulomb 

repulsion. From the positional and the orientational correlation functions, we have 
found an indication that the solid phase has a quasi-long-range (power-law corre- 
lated) positional order and a long-range orientational order. This implies that the 
long-range 1/r system shares the absence of the true long-range crystalline order 

qv ' at finite temperatures with short-range ones to which Mermin's theorem applies. 

We also discuss the existence of the "hexatic" phase predicted by the Kosterlitz- 
Thouless-Halperin-Nelson- Young theory. 
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1 Introduction 



More than 60 years ago, Wigner pointed out that an electron system will crys- 
tallize due to the Coulomb repulsion for low enough densities (Wigner crystal- 
lization) [1]. Although quantum effects play an essential role in a degenerate 
electron system, the concept of Wigner crystallization can be generalized to 
a classical case where the Fermi energy is much smaller than the thermal en- 
ergy. A classical two-dimensional (2D) electron system is wholly specified by 
the dimensionless coupling constant T, the ratio of the Coulomb energy to 
the kinetic energy. Here T = (e 2 /47rea)//c£T, where e is the charge of an elec- 
tron, e the dielectric constant of the substrate, a the mean distance between 
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electrons and T the temperature. For T ^C 1 the system will behave as a gas 
while for r > 1 as a solid. Experimentally, Grimes and Adams [2] succeeded 
in observing a transition from a liquid to a triangular lattice in a classical 2D 
electron system on a liquid-helium surface around T c = 137 ± 15, which is in 
good agreement with numerical simulations [4-8]. 

On the theoretical side, two conspicuous points have been known for 2D sys- 
tems: (i) Mermin's theorem dictates that no true long-range crystalline order 
is possible at finite T in the thermodynamic limit [9]. To be precise, the 1/r 
Coulomb interaction is too long ranged to apply Mermin's arguments directly. 
Although there have been some theoretical attempts [10,11] to extend the the- 
orem to the Coulomb case, no rigorous proof has been attained, (ii) A theory 
due to Kosterlitz, Thouless, Halperin, Nelson, and Young (KTHNY) predicts 
that a "hexatic" phase, characterized by a short-range positional order and a 
quasi-long-range orientational order, exists between a liquid phase and a solid 
phase [12]. Because the KTHNY theory is based on various assumptions and 
approximations, its validity should be tested by numerical methods such as 
a molecular dynamics (MD) simulation. Several authors have applied numer- 
ical methods to classical 2D electron systems, but they arrived at different 
conclusions on the KTHNY prediction [5-8]. 

In order to address both of the above problems, the most direct way is to 
calculate the positional and the orientational correlation functions, which is 
exactly the motivation of the present study [18]. 



2 Numerical Method 



A detailed description of the simulation is given elsewhere [17], so we only re- 
capitulate it. We consider a rectangular area with a rigid uniform neutralizing 
positive background in periodic boundary conditions. The aspect ratio of the 
rectangle is taken to be L y /L x = 2/V3, which can accommodate a perfect 
triangular lattice [3]. The Ewald summation method is used to take care of 
the long-range nature of the 1/r interaction. We have employed Nose-Hoover's 
canonical MD method [13,14] to incorporate temperature accurately. 

The system is cooled or heated across the transition with a simulated annealing 
method. The results presented here are for N = 900 electrons with MD runs 
with 30,000-110,000 time-steps for each value of T. The correlation functions 
and other quantities are calculated for the last ~20,000 time-steps of each run. 

Following Cha and Fertig [15], we define the positional and the orientational 
correlation functions from which we identify the ordering in each phase. The 



positional correlation function is defined by 
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where G is a reciprocal vector of the triangular lattice with the summation 
taken over the six G's that give the first peaks of the structure factor [inset 
(a) of Fig. 1]. The orientational correlation function is defined by 

J2 s ( r - l r i- r il)$s( r iMsfo). 
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where ip6( r ) — ~ Z)a' n ' e 6l6, " ( - r - ) with a (r) being the angle of the vector con- 
necting an electron at r and the a-th nearest neighbor with respect to a fixed 
axis. The summation is taken over n c nearest neighbors that are determined 
by the Voronoi construction [16]. 



3 Results and Discussions 



Let us first look at the positional and the orientational correlation functions in 
Fig. 1 for T = 200 and T = 160, for which the system is well in the solid phase. 
The positional correlation is seen to decay slowly, indicative of an algebraic 
(power-law) decay at large distances. The round-off in the correlation function 
around half of the system size should be an effect of the periodic boundary 
conditions. The algebraic decay of the positional correlation indicates that the 
2D electron solid has only a quasi-long-range positional order. Thus we have 
obtained a numerical indication that Mermin's theorem [9] remains applicable 
to the 1/r Coulomb interaction, which is consistent with the analytical but 
approximate results obtained in [10,11]. 

On the other hand, the orientational correlation rapidly approaches a constant, 
indicating a long-range orientational order. Therefore, while the 2D electron 
solid has no true long-range crystalline order, we can say that it has a topo- 
logical order. From a snapshot of the configuration [see inset (b) in Fig. 1], 
we can see that a long-range orientational order is preserved since defects 
(5- or 7-fold disclinations, etc) tend to appear as dislocation (5-7 combina- 
tion of disclinations) pairs, i.e., 5-7-5-7 disclination quartets that only disturb 
the orientational correlation locally. Here the coordination number is again 
determined from the Voronoi construction. 



Now we move on to the orientational correlation function (inset of Fig. 2) 
around the crystallization, which is obtained by cooling the system from a 
liquid to a solid. In between a short-range orientational order for T = 120 and 
a long-range one for T = 140, the correlation appears to decay algebraically 
at T = 130 with an exponent approximately equal to unity, which deviates 
from the upper bound of 1/4 predicted by KTHNY [19]. However, numerical 
difficulties arising from finite-size and finite-time effects prevent us from draw- 
ing any definite conclusion on the existence of the hexatic phase. Namely, we 
cannot rule out the possibility that the power-law behavior is an artifact of in- 
sufficient equilibration. In fact, the solid phase persists down to T = 130 when 
the system is heated from a solid, which is understandable if the solid-hexatic 
and the hexatic-liquid transitions are of first and second order, respectively. 

The KTHNY theory is based on a picture that the hexatic-liquid transition 
occurs through unbinding of disclination pairs. To examine if this is the case, 
we have calculated a defect-defect correlation function (Fig. 2), which we de- 
fine as a distribution of 7-fold coordinated electrons with respect to a 5-fold 
coordinated electron. The correlation function exhibits no qualitative differ- 
ence between T = 120 and T = 130, for which the disclinations are not tightly 
bound as compared with T = 140. If we look at a snapshot for T = 130 (Fig. 3), 
we see some domain structure as far as the present numerical simulation with 
finite-size and finite-time restrictions are concerned. A finite-size scaling will 
be an interesting future problem if the transition into the hexatic phase is of 
second order. 

We have also looked at a dynamical property, i.e., the power spectral density 
of the velocity (Fig. 4), which is related to the vibrational density of states and 
corresponds to the Fourier transform of the velocity autocorrelation function 
via Wiener-Khinchin's theorem. While the difference between the solid and 
the liquid phases appears in the spectrum around zero frequency, which is 
proportional to the diffusion constant (finite in the liquid or vanishingly small 
in the solid), we find a peak around the typical phonon frequency, which, 
curiously enough, persists even in the liquid phase. This indicates that the 
liquid has well-defined local configurations despite the short-range positional 
and orientational correlations. 
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Fig. 1. The positional (the upper frame) and the orientational (the lower frame) 
correlation functions for T = 200 and T = 160. The horizontal scale is the distance 
in units of the lattice constant of the triangular lattice. The structure factor for 
r = 200 is shown in inset (a), and a defect structure that is found for T = 160 
is shown in inset (b), where 5- fold (7- fold) coordinated electrons are marked with 
open (solid) circles. 



Fig. 2. Distribution of 7- fold coordinated electrons with respect to a 5- fold coor- 
dinated electron near the crystallization. The orientational correlation function is 
also shown in the inset. 



Fig. 3. A snapshot of the electron configuration for T = 130. 5-fold (7-fold) coordi- 
nated electrons are marked with open (solid) circles. 



Fig. 4. Power spectral density of the velocity for T = 120 and T = 140. The plot for 
r = 130 is similar to that for T = 120. The frequency scale is for n = 10 12 /m 2 . 
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